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The use of differential equations such as Eikonal, Hamilton-Jacobi and Poisson for the 
economical calculation of the nearest wall distance d, which is needed by some turbulence 
models, is explored. Modifications that could palliate some turbulence-modeling anomalies 
are also discussed. Economy is of especial value for deforming/adaptive grid problems. For 
these, ideally, d is repeatedly computed. It is shown that the Eikonal and Hamilton-Jacobi 
equations can be easy to implement when written in implicit (or iterated) advection and 
advection-diffusion equation analogous forms, respectively. These, like the Poisson 
Laplacian term, are commonly occurring in CFD solvers, allowing the re-use of efficient 
algorithms and code components. The use of the NASA CFL3D CFD program to solve the 
implicit Eikonal and Hamilton-Jacobi equations is explored. The re-formulated d equations 
are easy to implement, and are found to have robust convergence. For accurate Eikonal 
solutions, upwind metric differences are required. The Poisson approach is also found 
effective, and easiest to implement. Modified distances are not found to affect global outputs 
such as lift and drag significantly, at least in common situations such as airfoil flows. 
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Nomenclature 

turbulence modeling constant 

lift and drag coefficients 

wing chord 

nearest wall distance 

distance function 

array or grid point indices 

length scale function or modified Laplacian 

metric term in transformed equations or Mach number 

number of iterations 

exponent, weighting/biasing parameter or outward facing normal 

marker in marching front procedure for seed, trial or points to be later solved 
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residual for U equation 

Reynolds number based on wing chord 

front propagation velocity implied in Eikonal equation 

Cartesian velocity components for Eikonal front propagation 

Cartesian coordinates 

distance to nearest surface in ‘wall units’ 

grid expansion or under-relaxation parameter 

diffusion analogous coefficients in equations (3,6,7) 

constant in diffusion analogous coefficient evaluation (see Equation (10)) 

viscosity solution coefficient 

transformed coordinates 

dependent variable in differential wall distance or distance function equations 
domain for wall distance or distance function computation 
reduced domain for wall distance or distance function computation 
solid angle 


I. Introduction 

W all distances, d, are still a key parameter in many turbulence modeling approaches 12 and also in peripheral 
applications incorporating additional solution physics 3 4 . Also, their near wall iso-values can be used in grid 
generation 56 . Far field d contours can be used as a rapid means of evaluating computational interfaces on 
unstructured overset meshes having relative movements 7 . 

Surprisingly, with search procedures, the effort in calculating d can be significant. For example, even with Cray 
C90 class computers and time invariant meshes it can take 3 hours just to gain d in a large three-dimensional case 8 . 
Because of the expense and inter-block communication issues, in some codes the following dangerous 
approximations are made: 

i) computing distances down grid lines, not allowing for grid non-orthogonality; 

ii) computing d as the distance between a field point and the nearest wall grid point, instead of truly the surface; and 

iii) in multiblock grids, determining d on a purely block wise basis, ignoring the possibility that the nearest wall 
distance might be associated with another block. The latter can create large inaccuracies and also non-smooth, 
unhelpful-to-convergence, d distributions. In relation to point (iii), for overset grids the situation can arise where the 
same point in space has different equations depending on which block it is viewed from. In such a case, there is no 
reason why the solution should converge at all. 

Clearly, inexact, non-smooth or grid-blocking-dependent wall distances d will mostly be unhelpful. However, 
the deliberate modification of d to some distance function d can alleviate certain secondary turbulence model 

anomalies or extend modeling potential 9 . For example, d » d can alleviate the excessive influence of sharp convex 
features in the geometry on the turbulence model 14 . For comers or bodies/surfaces in close proximity, the increased 
multiple surface turbulence damping effect (see Mompean et al. 10 . Launder et al. 11 ) can be crudely modeled by 
setting d <d . 


A. Requirements of an ‘Ultimate’ Distance Function 

Following Ref. 11a preferable distance-function ( d ) behaviour is perhaps best captured in terms of the 
elemental solid angle, do). This is the angle subtended by a patch, on say a solid surface of radius of curvature R, a 
distance d from a field point. Hence Ref. 1 1 proposed 

i-Aj ,i) 

d n d 
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(here we correct the factor of 2 error in the equation of Ref. 1 1) where the continuous wall distance d can be found 
from say a nearest surface search (NSS). Here, the turbulence length scale function is L = 1, but for certain 
turbulence modeling requirements this may not be best 11 . Turbulence models are calibrated with a single flat surface 

mostly in mind with d « R or C0~ 2 k. In that case, from Equation (1), d = d. For several surfaces at about the 

same distance, the blocking effects are weakly additive and so d < d seems sensible. However, d should not greatly 
deviate from d. For instance, a reasonable effect, for two surfaces, would be to have the following harmonic related 
mean 1 Id" = 1 Id" +1 Id" where d t and d 2 are the distances to the two walls and n = 2 (in Ref. 11, n = 1 is 
suggested). For a channel, helpfully, considering the harmonic type mean with n = 2, CD = 4 k and for a flap cove 
CO > 2 7 T i .e. extra turbulence destruction is naturally introduced. 

The opposite situation is when the solid body is much smaller than its distance to the field point under 
consideration, or in other words the total solid angle it covers is much smaller than 1. An example is a thin wire, 
which clearly has a much weaker damping effect than a large flat surface at the same distance. We then need to have 
d » d. This also helps remedy the excessive modeled turbulence destruction that can be found around extreme 
convex features 3 4 . 

The above are all accuracy, or physics, considerations, and rather preliminary. Numerically, d field smoothness 
and computing speed are always desirable. This is more the key focus of the current paper. In addition, the accurate 
integration of Equation (1) with high-aspect-ratio grid cells near the wall is far from trivial, and the approach of 
directly solving (1) may be quite unattractive in practice for that simple reason. 


B. Differential d Equations 

Motivated by the expense of the d solution Sethian 5 considered the Eikonal equation below 


V0 = i+/i v> 


(2) 


seeking viscosity solutions where X — » 0. The dependent variable in equation (2), </), models propagating front first 
arrival times. The right hand side implies (away from shock-like features) the front has unit velocity, i.e. there is 
some velocity field with IUI = 1. This means first arrival times are equal to d. Fares and Schroder 1 essentially solve 

an Eikonal related equation for (jf x . To enhance modeling potential, the Laplacian has some control over d . The 
Eikonal equation with an explicit Laplacian, as below, is called a Hamilton-Jacobi (HJ) equation. 

N = i+r (3) 


Solving for <j>~' as in Ref. 1 could overcomplicate programming and requires an additional arbitrary length scale d 0 
to avoid infinite values at the wall. Also, in Ref. 1 no observation of the need for upwind differenced metric terms 
(see later) is made. The paper does not mention that the inverse d equation is connected to the Eikonal and HJ 
equations, and hence is amenable to specialized solution approaches. 

Spalding 4 approximately reconstructs d from solution of a more numerically benign Poisson equation 
(V 2 0 = - 1). The d(</>) reconstruction involves an auxiliary analytically derived equation. The Eikonal, unlike the 
Poisson approach is challenging to code 4 ' 5 . This is especially the case for unstructured grids. Hence its 
implementation in established industrial CFD solvers represents a significant code developer time investment and 
hence cost. Therefore, here, use of an Eikonal equation form, amenable to general geometry CFD code 
implementation is explored. The form is reminiscent of the Euler/Navier-Stokes equations (the key equations 
modeled in CFD solvers). Also, since the Poisson d method is easy to implement in industrial codes, this approach is 
further studied. For the current work, as the base CFD solver, the NASA CFL3D 12 program is primarily used. 

Operation counts show that differential approaches can be significantly faster than search procedures. For fixed 
mesh problems it is unusual for the search procedure to contribute significantly to the total solution cost. 
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Nonetheless, for increasingly-common non-stationary-mesh flow solutions the repeated search cost becomes highly 
significant. Then, differential approaches, which can make easy, safe use of previous d estimates, become especially 
attractive. The new d approaches presented here are best viewed in this moving mesh context. Although the 
approaches to be presented have successfully been tested for distorting meshes, only stationary mesh results are 
given here. Stationary meshes (where the starting/initialized d field is very different from the actual) most strongly 
test robustness and convergence traits and so seem appropriate. The equations considered are now more fully 
outlined. 


II. Implicit/Iterated d Equations 

In this paper, three different wall distance methods based on differential equations are used. These are: implicit 
Eikonal, implicit HJ and Poisson, defined below. The Eikonal (Equation (2)) is an exact d equation. When defining 
the vector 


U = W0 


(4) 


it can be rewritten in the following implicit (Eikonal!) advection analogous form 


U» W0 = \. 


(5) 


In the above, for convenience AV 2 0 where K — > 0 is omitted. The vector U corresponds to the front propagation 
velocity implied in the Eikonal equation where IUI = 1. With (4), the HJ (Equation (3)) can also be written in the 
following implicit (HJ!) form 


\]»V0 = \ + T(0)V 2 0 


(6) 


The positive function Y{0) is discussed later. Since (V 0\ = V • (0V 0)- 0V 2 0 and U = V0 the following conservative 
form of (6) is also possible 


v»(^u)=i+rv 2 ^ 


(7) 


where 


r = (r + 0 ) 


( 8 ) 


To gain Eikonali solutions, using an initial 0 field guess, equations (4) and (5) can be solved in an iterative sequence. HJ! 
solutions can also be made using the same approach. 

Along different lines, substitution of the Poisson based 0 distribution arising from solving Equation (6) with U = 0 
and T = 1 into 
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also gives distances 4 . This is called the Poisson method. The derivation of (9) assumes extensive (infinite) coordinates in 
the non-normal wall directions. Hence, unlike the Eikonal, d from Equation (9) is only accurate close to walls. However, 
turbulence models only need d accurate close to walls. 

A. Laplacian form and role 

Near a fine convex feature (for example a wire) for theoretical correctness accurate distances are needed so d = d . 
However, to prevent excessive far field influence d » d can be required further away from the feature 1 ' 4 . Adjacent to a 
convex feature TV 2 0»O. Therefore, the positive Laplacian inclusion in (6) has the desired effect of 

enlarging/exaggerating d(=0). Motivated by dimensional homogeneity, the need that as d — » 0 ,d = ribut V 2 </> — > °° 
suggests 

r = ed (io) 

where fis a constant. Clearly more ‘aggressive’ functions than (10) (e.g. T = e(-l + e d ) ) are possible but these are not 
explored. 

At concave comers V 2 <f> < 0 , hence cl « d. Therefore, with the Laplacian, the damping effects of ‘extra’ walls, 
discussed earlier, is naturally accounted for. Eikonal and hence Equation (5) distances will generally be discontinuous in 
gradient (this partly occurs where fronts effectively propagating from different walls collide). Therefore, the V 2 </> term 
smoothing in Equation (6) has the potential to enhance convergence. 

III. Numerical Modeling 

For brevity, the Equations (4-7) numerical modeling is described. Modeling for the much less challenging Poisson d 
method is discussed in Ref. 4. Equations (4-7) can be solved on curvilinear and unstructured grids. For curvilinear grids 
they must be transformed using the chain mle for differential calculus in say an £ rj system. When solving in this 
system, metric terms (M) such as £ r (=3£/3:t) , // , £ , and // are evaluated. This aspect will be discussed later. 

Equation (6-7) diffusive terms are discretized using second order central differences. The advection analogous 
derivatives in equations (5-7) are typically discretized, just considering an x coordinate direction for example, using the 
following 1 st order upwind type of approximation 



5 

American Institute of Aeronautics and Astronautics 



and h m = 1 for u. , > and m. , > 0 else « f l =0; and n. +1 = 1 for — m. +i > u. , and -M. +1 >0else n M = 0. The 
propagation velocities in the above are evaluated from (4). Where, for example 


</>, ~0,-l 

Av r ._, 


A T+1 


(13) 


Setting T = 0 in (7,8) shows that the conservation form of the Eikonal] equation still has a Laplacian. Here the view is 
taken that the discretized conservation form of the equation is being solved. However, to allow optimal Eikonali solution 
(using a hyperbolic advancing front approach described later), the extra Laplacian </IV 2 </> is ignored. The approximation 
of neglecting </>V 2 </> is equivalent to solving the non-conservation Equation (5) with the alternative approximation that 


U . = «. jM. 


+ n ,,u 


(14) 


The above equation corresponds to using offset difference-based velocity components. These correspond to those 
needed in the conservative Equation (7) form. Either viewpoint has no great accuracy implications. Where d needs to be 
accurate, for turbulence models, </iV 2 </> — >0. Numerical tests confirm the validity of the above arguments. 
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Figure 1: MF solution approach with one seed. 
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A. Metric Discretization 

The metric terms (M) must be carefully discretized. In CFL3D, for the i index direction M lA and M, +1 are first evaluated. 
The latter, for example, involves geometric data at the i and i+1 grid points. These M iA and M i+l difference based 
constructs are then, as with most CFD codes, averaged to give a single M; value to be used in the discretized transport 
equations. However, for strongly expanding near wall grid spacings this standard approach results in d overestimations 
for the Eikonali and HJi. To remedy this the following metric formulation is used 


M, = + n M M M 


(15) 


The standard CFL3D implementation corresponds to n = n M = 1/2 in the above. 

B. Stabilizing Measures 

To ensure stable solutions for the F.ikonah and HI) velocity clipping and diagonal dominance enhancement are tried. 
Both use the observation that, in 2D (two-dimensions), the exact U field should satisfy 


R = u 2 +v 2 -1 =0 

U I,J I,J 


(16) 


where u and v are r and y direction velocity components, respectively. Therefore, to improve diagonal dominance R u and 
R u </>i are added to both the discretized equation matrix diagonals and source temis, respectively. Based on Equation (16) 
the following velocity clipping is used 


u. j < 1, V, j < 1 


(17) 


The Eikonal equation does not permit a backwards-front movement 5 . The implication of this is that if in (6), for example, 
TV 2 0 < -1 a theoretical violation has occurred (the sign of the Equation (6) right hand side gives the front propagation 
direction). Hence, for HJj solutions the Equation (6) Laplacian (L) is modified to 


L = max - C, TV 2 ^] 


(18) 


where 0 < C < 1 . Here C = 1 is used. This is the theoretical correctness upper limit. Since anti-diffusion is associated 
with instability use of (18) should be viewed as a stability measure. As a further stability measure under-relaxation is 
used either through the following (/> = ( 1 - a)(f)° ,d + at/)"™' (where a is an under-relaxation factor and the ‘new’ and ‘old’ 
superscripts indicate iterative states) or use of a pseudo time temi. 

When solving the Eikonal equation explicitly, for each grid point, evaluation of the discretized matrix coefficients 
needs around 160 operations. The Eikonalj equation needs 90 operations but then iteration is also required (i.e. equations 
(4,5) must be contained in some iterative loop). Therefore, the key efficiency issue is the number of iterations N typically 
required. 


C. Simultaneous Equation Solution 

For Eikonali solutions either essentially crude global Gauss-Siedel (GGS) type iterations or a marching front (MF) 
approach are used. For HJi solutions the MF approach can be used to gain an initial d field and then GGS used. 
Alternatively, as for the Eikonalj, purely GGS based solutions can be made. To solve the Poisson equation essentially a 
GGS approach is used. 
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The MF approach essence is illustrated using Fig. 1. In this, a ‘wire’ is represented on a Cartesian grid using a single 
node point. Points surrounding this are considered to be either seed, trial or those to be later solved. These points are 
labeled through the marker variable ns with the following respective values 2, 1 or 0. To start, a seed (or multiple seeds) 
point is defined. The Eikonal; equations are applied N times (until convergence) to assign this seed point a cl value. The 
seed is represented in Frame (I) by the closed symbol. The Eikonalj equations are then solved for all immediate seed 
neighbors (trial points). These ns = 1 points are shown as open symbols. As shown in Frame (II), the trial point with the 
minimum cl is taken as the next seed. Now, for this point, ns = 2. Then, as shown in Frame (III), cl for all immediate 
neighbors to the new seed (for which ns ^ 2) are calculated. Hence the next seed point location is found. Frames (IV) and 
(V) show subsequent development stages for a clockwise moving circular front. In summary the procedure is as follows: 

a) Define a seed point; 

b) Iteratively solve for d at immediate neighbors to the seed (trial points) and 

c) Make the trial point with the minimum cl the next seed and return to (a). 

This procedure is continued until ns = 2 for all points in Q or a smaller domain/extent Q. . The latter domain is 
identified by the Fig. 1 dotted line. The Eikonal/MF method’s compatibility with smaller domains is useful. This is 
because, many turbulence models and DES (Detached Eddy Simulation) 2 only need d to a maximum of about l/3 rd 
the boundary layer thickness. In a practical system, all surface adjacent points can be taken as seeds and what is 
called an ‘active front’ produced. For best efficiency a heap-sort procedure is required. This has not been 
incorporated as part of this work. Cleary, based on operation counts, N < 2 is needed for the Eikonal, approach to 
match the Eikonal’s efficiency. However, a typical Eikonal I iV value is 4 (See Ref. 13) i.e. for stationary meshes the 
Eikonati equation has about twice the computational cost of the Eikonal. For moving meshes this difference greatly 
decreases. 


D. Initial/Starting Conditions 

For the Eikonal and III, 1 < </> < <=» in Q is found adequate. However, for MF solutions cp —> oo is used. Note, when 
using the MF approach to give a DES distance field cp = C DE sAj ( C DE s is a turbulence modeling constant) can be used 
where A, = max(Av,Ay,Az). The d computation will naturally stop when cl = C DES A,-. No DES d results are shown here. An 
example can be found in Ref. 13. For the Poisson approach <p~0 is adequate. 




Fig. 2: Badly distorted flat plate grid. Fig. 3: Predicted against actual datx- 0.5. 


E. Boundary Conditions 

Conditions on the domain boundaries are now described. At solid walls the following Dirichlet condition is applied 


<p = 0 (19) 
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At flow/far field boundaries 


^ = 0 


( 20 ) 


can be used, where n is the boundary nonnal co-ordinate. However, if Q is sufficiently large, (19) makes a 
stable, most computationally economical far field boundary condition. It is especially preferred for the Poisson 
method where it gives much faster convergence. 


IV. Results and Discussion 

The following geometries are considered: (a) Rat plate; (b) Single element airfoil (NACA4412); (c) Wing body and (d) 
Wing flap. Where d deviation (‘error’) values are given, these are relative to NSS distances. The deviation is defined 
using the equation below 


d „„ — d 
Error = 100(— ) 

d jvp. 


( 21 ) 


where d N ss is the distance coming from the nearest surface search. The d without the subscript in Equation (21) is the 
distance field for which the error is being evaluated. Row solutions use the Spalart and Allmaras 14 (SA) turbulence 
model. This probably has the most extreme d requirements of most standard turbulence models. It needs d accurate for 
around l/3 ld the boundary layer thickness. Hence it is a good candidate for d modeling tests. 




(a) x (b) 

Fig. 4: Distorted grid V^data: (a) IUI and (b) U vectors. 


A. Flat Plate (Case (a)) 

Initially, for testing just the Eikonali method, a flat plate is used. This is at y = 0, in a 2D square domain £2, with sides of 
unit length. First a grid (not shown) with strong grid expansion in the y direction is considered. Table 1 gives % d errors 
for different approximate geometric grid expansions a;, with offset and centered metrics. 
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OCy 

% d Error 

Offset 

metrics 

Centered 

metrics 

1.1 

0.48 

1.53 

1.7 

0.99 

8.25 

2.6 

3.15 

15.0 


Table 1. Eikonal; % d errors with 
offset and centered metrics. 


Clearly, (see Table 1) with centered differences and larger a, values serious errors arise. This is not surprising. Roache 15 
shows for accurate flow solutions ck,, < 1.3 is needed. The front propagation nature of the Eikonal; (and Eikonal) makes 
errors additive. Furthermore, centered differences are inconsistent with a hyperbolic propagating front problem. Hence, 
all the remaining Eikonal; results use offset metrics. 

Next, the effect of using a badly distorted grid is explored. Fig. 2 shows the badly distorted ‘curvilinear’ grid used in 
this study. This grid is intended to severely test robustness of the Eikonal; method. Fig. 3 plots predicted against actual 
distances (y) at .r=0.5 for the Fig. 2 grid. The full line gives the exact solution. The symbols give the distances predicted 
by the Eikonal; method. The figure shows that, even with an extremely distorted grid, reasonable distances can be 
obtained. 



Figure 5: Distorted grid plot of log (R„) 
against GGS iterations. 


Fig. 4 shows Fig. 2 grid velocity data. Frames (a) and (b) give IUI contours and U vectors, respectively. Considering the 
poor grid form the velocity magnitudes and direction are surprisingly good. Fig. 5 gives a Fig. 2 grid plot of log (R u ) 
against GGS iterations. Through its monotonic nature the plot reflects the stable convergence. 


B. Single element airfoil (Case (b)) 

The Poisson, Eikonal; and NSS procedures are now compared for a NACA4412 airfoil. The attack angle is 13.87°, 
M = 0.2 and Re = 1.52 x 10 6 (based on the wing chord, c). Fig. 6 shows the 3 zone circa 25,000 cell overset grid used as 
part of the studies. Fig. 7 gives Eikonal; U vectors. Correctly, the vectors are surface normal. Their magnitude is close to 
unity. Overset Poisson solutions are made for Q and also just the inner zone region ( L2 ) labeled in Fig. 6. Obviously, 
the latter solution is most computationally efficient. At the Q. boundaries the Dirichlet condition 0 = 0 is used. Fig. 8 
shows d contours representing in the following respective frames the: (a) NSS; (b) Eikonal;; (c) Poisson and (d) 
Q Poisson results. 
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Fig. 9 gives v + <400, d deviation histograms. Frame (a) is for the Eikonab equation. The average deviation is 0.5 %. 
Frame (b) is the overset Poisson Q. domain solution. Although not shown by the histograms, the average d deviation for 
the Poisson is significantly higher (2.5%) than that for the Eikonal equation. 
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Fig. 10 gives a zoomed in view of trailing edge region Poisson cl deviation contours. As can be seen, a key d 
overestimation zone is the sharp convex trailing edge geometry region. This overestimation zone, which will also be 
produced with the H.^ equation, is potentially desirable (see later). 



C L 

c D 

Eikonal[ 

1.704 

0.03466 

NSS 

1.698 

0.03496 

Poisson 

1.713 

0.03543 


Table 2. C L and C D ’s for different d fields. 

Table 2 gives lift ( C L ) and drag (C D ) coefficients. Some minor differences are evident. For surface pressure plots and 
boundary layer velocity profiles lines/results effectively overlay. Fig. 11 gives SA, trailing edge region, turbulent 
viscosity (//,) contours for the NSS (Frame (a)), and Poisson (Frame (b)) approaches. The larger x/c > 1 Poisson 
distances, reduce the modeled turbulence destruction thus slightly increasing //, by just under 5%. Tests for flow over a 
thin wire 4 suggest turbulence levels for the exaggerated distances provided by the HJ and Poisson methods, are more 
realistic than those provided by the NSS. 

Clearly, for the current case, turbulence intensities downstream of the airfoil have insignificant influence on the 
parameters of interest in a design context. However, for multiple element airfoils it is not inconceivable that more 
significant solution differences could arise. For example, these can be caused by a peak in modeled turbulence energy 
convecting close to the center of the leading edge of a downstream element. Then small changes in the turbulence energy 
peak’s position can give rise to different solutions. 

The Poisson and NSS methods were also used as part of a single airfoil ‘aeroelasticity’ analysis comparing with 
AGARD Case 3 measurements 16 . The Poisson equation shows fast convergence and gives instantaneous C L , C D and 
moment coefficient values very similar to NSS. These results are reported in Ref. 13. 






(d) 


Figure 8: Case (b) (NACA4412) d contours: (a) NSS; (b) 
Eikonah; (c) Poisson and (d) inner block Poisson. 


12 

American Institute of Aeronautics and Astronautics 






C. Wing body (Case (c)) 

For this case the angle of attack is 2.87°, M = 0.802 and Re = 13.1 x 10 6 (based on c). The single block grid, for which 
the body surface zone is shown in Fig. 12, has around 0.9 million cells. For this case the Poisson, Eikonalj and NSS 
approaches are tested. Fig. 13 gives the y + < 400-region d deviation histogram for the Eikonal! and Poisson methods. 
The average ‘errors’ are 3.13 % and 3.14 %, respectively. The Poisson’s tendency to over predict d is evident in the 
histograms. Table 3 gives Cl and Cd values for the different d fields. For this more complex geometry case, the average 
differences between C L and Cd is lower than for Case (b). 


D. Wing Flap (Case (d)) 

For a multi-element airfoil case, the wing and flap angles are 2° and 40°, respectively. The wing-flap gap is 0.6% of c. 
The Reynolds number is 23 x 10 6 , based on c, and M = 0.18. A small part of the highly stretched 10 block, 
approximately 0.9 million-cell grid is shown in Fig. 14a. For this case the Poisson, NSS and Eikonal! methods are tested. 
Fig. 14b gives Eikonalj U vectors. Fig. 15 gives a zoomed in view of area A. Correctly, the vectors are surface 
normal. Magnitudes are close to unity. Fig. 14c gives Eikonali d contours. Since they are virtually identical, NSS 
contours are not shown. Poisson contours can be found in Ref. 4. 
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Figure 9: Case (b)y+ < 400 d deviation histogram: 
(a) Eikonal! and (b) Poisson. 
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Figure 10: Case (b) (NACA4412) Poisson d deviation contours. 
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Figure 11: Turbulent viscosity contours: (a) NSS and (b) Poisson. 



C L 

c D 

Eikonal, 

0.6633 

0.04839 

NSS 

0.6632 

0.04855 

Poisson 

0.6635 

0.04842 


Table 3. Case (c) C L and C B ’s 
for different d fields. 



Figure 12. Case (c) (wing body) surface grid. 
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Fig. 16 gives y + < 400 cl deviation histograms. Frame (a) is for the Eikonal;. Frame (b) is for the Poisson. The 
average Frame (a) deviation is 2.1 %. The Poisson equation has a 1% deviation. This low Poisson ‘error’, relative to 
the Eikonal; equation seems to contradict the histogram evidence. The contradiction is because the Eikonal equation 
has some large deviations at just a very few points. Clearly, the Frame (a) Eikonal deviation distribution has the 
least spread. For the Poisson method (see Frame (b)), it can be seen there is a slight tendency to over predict d. As 
can be inferred from Table 4, Poisson method lift and drag coefficients are within 0.05% of those for the NSS 
procedure. For the Eikonal, equation the deviation is a little greater but not that significant. 

Pleasingly, complex geometry tests show 1 1 J , GGS convergence rate is similar to the Eikonal;. Hence the 1 1 J, 
Faplacian does not seem that important for stability, its value being more in offering the capability of controlling 
wall distances near convex or concave surfaces. 



c L 

c D 

Eikonal; 

2.815 

0.0582 

NSS 

2.810 

0.0584 

Poisson 

2.811 

0.0584 


Table 4. Case (d) C L and 
Cfl’s for different d fields. 



% Error 



Figure 13: Case (c) (wing body) y+ < 400 d deviation 
histogram: (a) Eikonal; and (b) Poisson. 
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Figure 14: Wing flap geometry near surface information: 
(a) grid; (b) U vectors and (c) d contours. 



Figure 15: gives a zoomed in view of area A. 
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d over prediction 


d under prediction 
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- 10.0 - 5.0 0.0 5 0 10.0 

% Error 



% Error 


Figure 16: Case (d) (wing flap) y+ < 400 d deviation 
histogram: (a) Eikonali and (b) Poisson. 


Y. Conclusions 

Implicit forms of the Eikonal (Eikonal!) and Hamilton-Jacobi (HJj) wall distance equations were presented. These 
are reminiscent of advection and advection-diffusion equations. Because of this, the Eikonali and 1 1 J| were found to 
be relatively easy to implement in an established industrial CFD solver. The implicit d equations were found to have 
robust convergence. Geometries studied included single and two element airfoils and a wing-body configuration. 
For [I.I|/Eikonal[ accuracy, offset metric differences are required. A simple, approximate, Poisson-equation-based 
distance approach was also found effective. Since it did not require offset metric evaluations it was easiest to 
implement. The sensitivity of flow solutions to plausible modifications of d was explored. Results were not greatly 
affected by wall distance traits, partly because peculiar geometries such as thin wires were not involved. 
Downstream of the airfoil trailing edge, intensities were altered by just under 5 %. 
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